// Trabalho de MEF 
// Problema bidimensional com calor sobre uma peça metálica
// 27/06/2012

clear;
clc;

// -------------------- Funções -------------------- //

// largura, altura, número de elementos na horizontal, na vertical, ponto não desejado
function [LG,EQ,q,vq,XY,npre] = montaMalha(w, h, numx, numy, p, v0)

    // Matriz LG nxm dos nós locais para cada elemento finito.
    // Sendo n o número de elementos e m o número de nós em cada elemento.
    // LG = [[1, 2, 7, 6]; [2, 3, 8, 7]; [3, 4, 9, 8]; [4, 5, 10, 9]];
    LG = zeros(numx*numy-4,4); // Os quatro elementos ao redor de p não entram
    
    // Vetor de correspondência dos elementos finitos para a matriz global.
    // Valores zerados significam que o elemento daquela posição no vetor 
    // não contribuirá na matriz pois possui valor prescrito.
    // EQ = [1, 2, 3, 4, 0, 5, 6, 7, 8, 0];
    EQ = zeros(numpx*numpy-1,1); // Não se inclui o p
    
    // Matriz coordenadas dos nós globais.
    // XY = [[0,0]; [1,0]; [2,0]; [3,0]; [4,0]; [0,1]; [1,1]; [2,1]; [3,1]; [4,1]];
    XY = zeros(numpx*numpy-1,2); // idem...
    
    // Vetor correspondente aos valores prescritos de cada nó.
    // q = [0, 0, 0, 0, val_pre, 0, 0, 0, 0, val_pre];
    q = zeros(numpx*numpy-1,1); // idem...
    vq = zeros(numpx*numpy-1,1); // idem... - máscara para o vetor q
    
    // Tamanho do passo
    hx = w / numx;
    hy = h / numy;
    
    numpx = numx+1;
    numpy = numy+1;
    
    metadex = -w/2;
    metadey = -h/2;
    
    // Preenchendo o vetor q com os valores prescritos
    q(p-numpx-1) = v0;
    q(p-numpx) = v0;
    q(p-numpx+1) = v0;
    q(p-1) = v0;
    q(p) = v0;
    q(p+numpx-2) = v0;
    q(p+numpx-1) = v0;
    q(p+numpx) = v0;
    npre = 8;
    
    vq(p-numpx-1) = 1;
    vq(p-numpx) = 1;
    vq(p-numpx+1) = 1;
    vq(p-1) = 1;
    vq(p) = 1;
    vq(p+numpx-2) = 1;
    vq(p+numpx-1) = 1;
    vq(p+numpx) = 1;
    
    // Número de pontos = número de elementos + 1
    ind = 1;
    for i = 0 : numy-1
        for j = 0 : numx-1
            p1 = i*numpx+j+1;
            p2 = p1+1;
            p3 = (i+1)*numpx+j+2;
            p4 = p3-1;
            
            if p1 <> p & p2 <> p & p3 <> p & p4 <> p then
                
                if p1 > p then
                    p1 = p1-1
                end
                
                if p2 > p then
                    p2 = p2-1
                end
        
                if p3 > p then
                    p3 = p3-1
                end
                
                if p4 > p then
                    p4 = p4-1
                end
                
                LG(ind,:) = [p1, p2, p3, p4];
                ind = ind + 1;
            end
            
        end
    end
    
    ind = 1;
    for i = 0 : numpy-1
        for j = 0 : numpx-1
            if i*numpx+j <> p then // Não queremos o ponto p, que corresponde à carga
                XY(ind,:) = [j*hx, i*hy];
                
                if XY(ind,1) == 0 | XY(ind,1) == h | XY(ind,2) == 0 | XY(ind,2) == w  then
                    vq(ind) = 1;
                    npre = npre + 1;
                end
                ind = ind + 1;
            end
        end
    end
    
    ind = 1
    for i = 1 : length(vq)
        if vq(i) == 0 then
           EQ(i) = ind;
           ind = ind+1;
        end
    end
endfunction

// --------------------------------------------------
// Função auxiliar para mostrar a situação da malha 
// em relação aos pontos prescritos e com a solução 
// em cada ponto (nó) da malha.
// --------------------------------------------------

function [] = mostraMalha(var, numx, numy, p, tipo, solucao)
    
    if tipo == "pre" then
        ind = 1;
        printf("\n");
        for i = 1 : length(var)
            if i == p then
                printf("- ");
                ind = ind + 1;
            end
            
            printf("%d ", var(i));
            
            if modulo(ind,numx+1) == 0 then
                printf("\n");
            end
            ind = ind + 1;
        end
        printf("\n");
    elseif tipo == "sol" then
        
        indx = 1;
        indy = 1;
        indvq = 1;
        indsol = 1;
        
        for i = 0 : numy
            indx = 1;
            for j = 0 : numx
                if i*numpx+j+1 == p then
                    var(indy, indx) = -1;
                    indvq = indvq - 1;
                elseif vq(indvq) == 1 then
                    var(indy, indx)= q(indvq);
                else
                    var(indy, indx) = solucao(indsol);
                    indsol = indsol + 1;
                end
                indx = indx + 1;
                indvq = indvq + 1;
            end
            indy = indy + 1;
        end
        
        disp(var);
    end
endfunction

function [saida] = phi(i,xi,eta)
    if i == 1 then
        saida = (1 - xi) * (1 - eta) / 4;
    end
    
    if i == 2 then
        saida = (1 + xi) * (1 - eta) / 4;
    end
    
    if i == 3 then
        saida = (1 + xi) * (1 + eta) / 4;
    end
 
    if i == 4 then 
        saida = (1 - xi) * (1 + eta) / 4;
    end
endfunction

function [saida] = dphi(i,var,xi,eta)
    if i == 1 then
        if var == "xi" then
            saida = -(1 - eta) / 4;
        end
        if var == "eta" then
            saida = -(1 - xi) / 4;
        end
    end
    if i == 2 then
        if var == "xi" then
            saida = (1 - eta) / 4;
        end
        if var == "eta" then
            saida = -(1 + xi) / 4;
        end
    end
    if i == 3 then
        if var == "xi" then
            saida = (1 + eta) / 4;
        end
        if var == "eta" then
            saida = (1 + xi) / 4;
        end
    end
    if i == 4 then
        if var == "xi" then
            saida = -(1 + eta) / 4;
        end
        if var == "eta" then
            saida = (1 - xi) / 4;
        end
    end
endfunction

// --------------------------------------------------
// Função jacobiano que utiliza a matriz de derivadas gerada pela função N
// a fim de fazer a correspondência de um elemento finito genérico para um 
// elemento que permita um cálculo de maior facilidade.
// --------------------------------------------------

function [invJ, detJ] = jacobiano(e, xi, eta)
    nos = LG(e,:);
    
    for i = 1 : nnos_el
        matriz_coord(i,:) = XY(nos(i),:)    
    end
    
    J = N(xi, eta) * matriz_coord;
    
    invJ = inv(J');
    detJ = det(J);
endfunction

function [saida] = N(xi, eta)
    N = zeros(2, nnos_el);
    
    for i = 1 : nnos_el
        N(1,i) = dphi(i, "xi", xi, eta);
        N(2,i) = dphi(i, "eta", xi, eta);
    end
    
    saida = N;
endfunction

// ---------- Funções K e F local ---------- //

function [saida] = montaKLocal(e)
    M = zeros(4,4);
    Ke = zeros(4,4);
    
    for i = 1 : num_pg
        [G,detJ] = jacobiano(e, PG(i,1), PG(i,2));
        _N = (N(PG(i,1), PG(i,2)));
        M = _N' * G' * Q * G * _N * detJ .* W(i);
        
        Ke = Ke + M;
    end
    
    saida = Ke;
endfunction

// Em nossa modelagem, não há fonte.
function [saida] = montaFLocal(e, k)
    nos = LG(e,:);
    parc = zeros(nnos_el, 1);
    
    for i = 1 : nnos_el
        matriz_coord(i,:) = XY(nos(i),:);
    end
    
    for a = 1 : nnos_el
        for b = 1 : nnos_el
            parc(a) = parc(a) - k(a,b) * q(nos(b));
        end
    end
    
    saida = parc;
endfunction

// --------------------------------------------------

// Monta matriz K global e vetor F global.
function [saida,saida2] = montaGlobals()

// Cria a matriz global.
KGlobal = zeros(nnos-npre, nnos-npre);
FGlobal = zeros(nnos-npre,1);

// Percorre cada elemento da superfície e gera sua matriz local.
for e = 1 : nel
    Klocal = montaKLocal(e);
    //disp(Klocal);
    Flocal = montaFLocal(e, Klocal);
    //disp(Flocal);
    for a = 1 : nnos_el
        for b = 1 : nnos_el
            i = EQ(LG(e,a));
            j = EQ(LG(e,b));
            if (i <> 0 & j <> 0) then
                KGlobal(i,j) = KGlobal(i,j) + Klocal(a,b);
            end
        end
        if (i <> 0) then
            FGlobal(i,1) = FGlobal(i,1) + Flocal(a);
        end
    end
end

saida = KGlobal;
saida2 = FGlobal;
endfunction

//-----------------------------------------
//---------- Início do programa ----------
//-----------------------------------------

// Coordenadas da malha
w = 10;
h = 10;

// Número de elementos
numx = 10;
numy = 10;

// Número de pontos (nós)
numpx = numx + 1;
numpy = numy + 1;

// Ponto não desejado (carga)
p = 61;

// Valor prescrito para a superfície equipotencial
v0 = 10;

// npre = número de nós com valores prescritos
[LG,EQ,q,vq,XY,npre] = montaMalha(w,h,numx,numy,p,v0)

nnos = length(EQ); // Número de nós
nel = length(LG(:,1)); // Número de elementos
nnos_el = length(LG(1,:)); // Número de nós por elemento

num_pg = 4; // Número de pontos de Gauss
c_pg = sqrt(3)/3;

// Matriz com os pontos de Gauss
PG = [[-c_pg,-c_pg]; [c_pg,-c_pg]; [c_pg,c_pg]; [-c_pg,c_pg]];

W = [1,1,1,1]; // Peso dado a cada ponto de Gauss

// Matriz de características do material
Q = [[1,0]; [0,1]];

printf("\nMatriz dos pontos prescritos:\n");
mostraMalha(vq, numx, numy, p, "pre");
//printf("Matriz EQ:\n");
//mostraMalha(EQ, numx, numy, p, "pre");

[K,F] = montaGlobals();
//printf("\nMatriz K:");
//disp(K);
//printf("\nVetor F:");
//disp(F);

solucao = K \ F;

printf("\nSolucao:");
//disp(solucao);
malhaRes = zeros(numx+1,numy+1);
mostraMalha(malhaRes, numx, numy, p, "sol", solucao);
